Deformable self-propelled particles 
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A theory of self-propelled particles is developed in two dimensions assuming that the particles 
can be deformed from a circular shape when the propagating velocity is increased. A coupled set of 
equations in terms of the velocity and a tensor variable to represent the deformation is introduced 
to show that there is a bifurcation from a straight motion to a circular motion of a single particle. 
Dynamics of assembly of the particles is studied numerically where there is a global interaction such 
that the particles tend to cause an orientational order. 
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Individual and collective dynamics of self-propelled 
objects are one of the fundamental problems in non- 
equilibrium physics. There have been a number of stud- 
ies for many years based on both deterministic dynamics 
0, and stochastic dynamics [1, 0, H, H, @, H|. Re- 
cently the hydrodynamics effects with and without in- 
ternal degrees of freedom have been investigated inten- 
sively [H, [l(| EH Formulation of a droplet or a 
vesicle moving under the interaction with substrate has 
been studied 13 L 3, [HI]. Experiments of reaction-driven 
propulsion llrl Il7ll as well as theory [l8[ and computer 
simulations [19, [2Q] have been carried out recently. Ex- 
citable reaction diffusion systems have also exhibited self- 
propelled domai ns J2ll| and have been analyzed theo- 
retically [H, [13, l24|. These studies are motivated, on 
one hand, to understand macroscopic self-organization 
far from equilibrium and, on the other hand, to clarify 
the mechanism of molecular machines in mesoscopic or 
nanoscopic length scales. 

All of the above studies except for the model with in- 
ternal degrees of freedom such as a dimmer model [HJ] 
assume that the particle shape is unchanged during the 
motion. In reality, however, many self-propelled objects 
may change the shape depending on the velocity or the 
interaction with other objects. Biological systems such as 
living cells are one example [26[. The excitable reaction 
diffusion system provides another example. The numeri- 
cal results in ref. [21] clearly show that a two-dimensional 
excited pulse changes its shape as the velocity increases. 

The purpose of the present Letter is to investigate the 
individual and the collective motions of self-propelled de- 
formable particles. It will be shown that a competition 
between a circular motion of individual particles and an 
orientational order of their shape due to a global inter- 
action causes a rich variety of collective dynamics. 

In order to represent the motion of a two-dimensional 
deformable particle with a coupling between the velocity 
of the center of the gravity and the deformation from a 
circular shape, we introduce two sets of variables. One is 



the velocity v = U2) of the center of gravity and the 
other is a tensor S a p with (a,/? = 1,2). Suppose that 
the deformation is weak and each particle deforms into 
an elliptical shape. We define a unit vector n along the 
long axis as shown in Fig. [T] The tensor S a p which is 
the same as the nematic order parameter [27| is given in 
terms of the components of n by 



S a p = s [ n a np - -S a p 



(1) 



which is normalized as TrS* = and the positive constant 
s represents the degree of deformation from a circular 
shape. 
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FIG. 1: Self-propelled particle with the velocity v and the 
unit normal n along the long axis. 



The time-evolution equations for v and S a p can be 
derived by a symmetry argument as follows. First of all, 
the velocity should obey in its simplest case 



dt 



jv a - \v 2 \v a - aS a pvp 



(2) 



where the repeated indices imply summation. The con- 
stant 7 is assume to be positive. The first and the sec- 
ond terms are the same as those of an active Brownian 
particle 0, SB]- The combination S a pvp with a scalar 
constant a as the last term of © is the simplest way to 
make a vector variable in terms of v and S a p. Similarly, 
by symmetry, a tensor v a vp — \\v \8 a p should enter into 
the time-evolution equation for S a p 



dt 



S a p = -nS a p + b(v a vp~^\v 2 \S a pj (3) 
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where n > and b are constant. Note that ([3|) satisfies 
TrS = 0. 
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By writing the velocity as v\ = vcos<ft, vi = usin^> 
and the unit normal as n\ — cos 9 and n 2 = sin#, eqs. 
(J2j) and ([3]) are rewritten as 

^-v = i! (7 - v 2 ) - ^asv cos 2(6* - tf>) (4) 



d_ 

db 
d_ 
dt 



-ft = —-as sin 2(9 — (/)) 

s = —Ks + v 2 bcos2(9 
v 2 b 



2s 



■sin2( 



(5) 
(6) 
(7) 



where v and s should be positive. From eqs. ([5]) and lj7]). 
we obtain for ip = 9 — <j) 

= _-(-a S + -)sm2^ (8) 

It is noted that only the difference V' of the two angles en- 
ters into the time-evolution equations due to the spatial 
isotropy. 

The above set of equations has a stationary solution. 
Without loss of generality, we may assume that the ve- 
locity is directed along the x-axis, i.e., <fi = 0. From eqs. 
([7]) and (JSJ, one notes that the angle 9 is given either by 
9 = or by 9 — ir/2 depending on the sign of b. That 
is, when b is positive, 9 — 0, and the steady value of the 
amplitudes is given by v 2 — 7/(1 + B) and sq = bv^/n 
where B = ab/(2n). Therefore the long axis of the ellip- 
tical particle is parallel to the propagating direction. On 
the other hand, when b is negative, the solution is given 
by 9 = it/2, v 2 = 7/(1 + B) and s = —bv 2 /n. In this 
case, the long axis is perpendicular to the velocity vector. 
In order to avoid a singular behavior of v Q for B < — 1, 
we hereafter assume ab > (B > 0). 

The linear stability of the stationary solution can be 
examined by a standard method. It is readily shown that 
the linearized equations for the amplitudes ^ and © 
are decoupled with those of the angle variables and any 
instability does not appear from the former. It turns out 
that the stability is solely determined by eq. (|8]). The 
stationary solution is stable when the coefficient of rhs 
of eq. ([5]) is negative. This gives us the condition that 
abv 2 , < k 2 and hence the threshold that 



7 = 7c = ^ + 2 



(9) 



When 7 > 7c, the stationary straight motion becomes 
unstable. Note that this is valid for both 9 = and 9 = 
7r/2. This bifurcation threshold is indicated on the 7— k 
plane in Fig. [5J The instability occurs when the velocity 
exceeds the critical value defined by v c = K/(ab) 1 ^ 2 . We 
assume that the value of v c is sufficiently small. This 
is required because the deformation is expected to be 
enhanced for larger velocity whereas the deviation from 
the circular shape in our model is expressed only by S a p 
with the mode 29 ignoring the higher harmonics. 
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FIG. 2: (a) Stability diagram on the 7 — k plane for ab = 0.5. 
(b) The velocity v (full line) and the frequency ui (broken line) 
as a function of 7 for a = —1.0,6 = —0.5 and n = 0.5. The 
bifurcation occurs at 7 C = 0.75. 



Now we investigate the dynamics of a self-propelled 
particle when the steady solution vq and sq of a straight 
motion becomes linearly unstable. Here we make an 
ansatz that there occurs a circular motion and explore 
the stability both numerically and analytically. To this 
end, we put v = v r , s = s r and 9 = wt + Q/2 and 4> = cut. 
Substituting these into eqs. ^ , (O , © and (O, one 
obtains after some algebra 



7-5 



(10) 



and s 2 = bv 2 /a, cos( = n/as r and lo 2 = (ab/4)(v 2 — v 2 ). 
Equation (jlpp requires that 7 > k/2 but it is found from 
eq. §§§ that this is automatically satisfied when 7 > 
7 C . It should also be noted that this inequality is the 
condition for the existence of the phase C- The frequency 
lo continuously increases from zero at v r = v c . The radius 
of the circular orbit is given by R = v r /uj. Comparing 
the stationary solution Vq = 7/(1 + B) with v r given by 
eq. (fTO)) , one notes that v is continuous at the bifurcation 
point 7 = 7c- The 7— dependence of v and lo is shown in 
Fig. db). 

In order to study the dynamics near the bifurcation 
7 = 7c in more detail, we make a reduction of the vari- 
ables. It is noted that the coefficient of eq. ([S]) vanishes 
at the bifurcation point. Therefore, the variable ip is slow 
near the stability threshold compared to other variables 
v and s. This allows us to eliminate these variables by 
putting dv/dt = ds/dt = in eqs. ^ and © so that eq. 
(IH1) becomes 



dV 



F(i>) 



7 - 7c + (k/2 - 7)(1 - cos 2 2ip) 



-27c + k(1 - cos 2 2i(j) 



tan 2i/j 
(11) 



The function F{ip) is an odd function since both clock- 
wise and counter clock-wise circular motions are equally 
possible. It is readily verified that eq. (jTTJ) for — it/ 2 < 
tjj < it 1 2 is monostable for 7 < 7 C whereas it becomes 
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bistable for 7 > j c . Therefore this is a pitchfork bifurca- 
tion. 
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FIG. 3: Trajectory of a particle for (a) K = K/k = 0.14 
and for the time interval St = 5000, (b) K — 0.16 and for 
St = 5000, (c) K = 0.28 and for St = 5000, and (d) K = 0.32 
and for St > 2000. Note that (c) and (d) cover wider area 
compared with (a) and (b). 



The linear stability of the circular motion is studied nu- 
merically by evaluating the eigenvalues of the linearized 
equations of ||3J), ^ and ([5]). It turns out that the sta- 
bility threshold coincides with the stability threshold of 
the straight motion given in Fig. [2] This implies that 
in the region where the straight steady motion is stable, 
the circular motion is unstable and vice verse. 

It is mentioned here that a model of self-propelled par- 
ticle which exhibits a circular motion has been introduced 
quite recently by Teeffelen and Loewen [29(. However 
their model does not have a bifurcation from a straight 
motion to a circular motion shown above. 

Now we explore the collective motion of the deformable 
particles. We introduce a kind of a global coupling such 
that 



d o(«) 
dt af3 



-kS 



(n) 
aB 



b [ v^v[ n) 



(12) 



where S = (^/N)J2 n ^af3 an ^ K ^ s the coupling con- 
stant. The superscript n indicates the n-th particle. The 
total number of the particles is denoted by N. Equa- 
tion (fl"2"f should be coupled with eq. ((2]) in which v a and 



S a n are replaced, respectively, by S^I and S y ' l J. The 
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An) 



'a/3 ■ 



last term in eq. (fT2"|) with a positive value of K implies 
that the ellipsoidal particles tend to make an orienta- 
tional oder. It is mentioned that an oricntational order 
has been studied in self-locomoting rods which are not 
deformable [H. 




(b) 



40 



20 







xxxxxxxxxxscxXx*?' 



0.1 0.15 0.2 0.25 0.3 0.35 

K 



0.1 0.15 0.2 0.25 0.3 0.35 

K 



FIG. 4: (a) Lyapunov exponent as a function of the inter- 
action strength K. (b) Diffusion constant D (plus) and the 
coefficient B (cross) in eq. (|15[) as a function of the interac- 
tion strength K. These are obtained from 100 independent 
runs. 



We have investigated the motion of the particles nu- 
merically in the situation such that individual particle 
undergoes a circular motion when K = 0. The parame- 
ters are fixed as a = —1.0, b = —0.5, n = 0.5 and 7 = 1.0 
whereas the interaction strength K = K/k is varied. The 
number of the particles is chosen as N = 30 in the most 
of numerical simulations. 

Figure [3] displays the trajectory of a particle for four 
different values of K. When K is small as K = 0.14, the 
motion is almost circular and localized as in Fig. EJa). 
However a transition to a delocalized state occurs around 
at K — 0.15. A random drift motion appears as in Fig. 
[3th). For larger values of K, a random but ballistic mo- 
tion becomes dominant as shown in Fig. [21c). When K 
exceeds 0.30, each particle moves straightly as displayed 
in Fig. 05d) . The direction of the velocity of each particle 
is completely random. 

In order to analyze the above behavior of the collec- 
tive motion, we have evaluated the Lyapunov exponent 
L defined through the relation 

N 

exp(L(i-t )/r) = £ - v^ n \t)f 

n=l 

-&ct{S {n) {t)-S {n) {t))] (13) 

where r = 2-k/lj — %^/2tt w 35.5 for the present set of 
the parameter values. Small perturbations of the order 
of 10~ 4 are introduced for the hat-marked variables at 
t = 3200. The exponent L is obtained by changing the 
value of K. The results are summarized in Fig. [4{a). 
The Lyapunov exponent is almost zero for K < 0.14 
and K > 0.3 whereas it is definitely positive for 0.15 < 
K < 0.29. Therefore we conclude that the motion is 
chaotic in this parameter regime. In Fig. HJa), the data 
of 100 independent runs are plotted. The vertical bars 
indicate the typical scatter of the data. The reason why 
the exponent is not negative but zero in the non-chaotic 
region is that there is a zero eigen-mode in the time- 
evolution equations due to the rotational invariance. 
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In order to quantify the distinct dynamics in these pa- 
rameter regions, we have calculated the averaged mean 
square displacement for all 30 particles defined by 



W 



N 



1 

N ^ 



,(») 



(t')dt' 



1/2 



(14) 



Since there is a crossover from a random drift to a ballis- 
tic motion of a particle as is evident in Fig. [3j we assume 
the time-dependence of W for sufficiently large values of 
t as 



W = 2D(t/r) 1/2 + Bt/r 



(15) 



The diffusion constant D defined through the relation 
(fT5)l is plotted in Fig. 0Jb) as a function of the interaction 
strength K. It is clearly seen that the diffusion constant 
starts to increase in the chaotic region and becomes zero 
for K > 0.3 whereas the ballistic motion becomes domi- 
nant and the value of B increases for K > 0.25 as shown 
in Fig. [U(b). The meaning of the diffusion constant is 
less appropriate for these large values of K. 

In summary we have introduced and studied a kinetic 
model for deformable self-propelled particles. It is found 
that an isolated particle exhibits a bifurcation such that a 
straight motion becomes unstable and a circular motion 
appears. It is noted that this is not an ordinary Hopf 
bifurcation and that the frequency of the circular motion 
increases continuously from zero at the bifurcation point. 
As indicated by eq. (fTTj) . the bifurcation is a pitchfork 
type accompanied with a zero mode due to the spatial 
isotropy. 



Assembly of the deformable particles exhibits complex 
dynamics where a global interaction is introduced in such 
a way that the particles deformed in an elliptical shape 
tend to have an orientational order. In the situation that 
each particle undergoes the circular motion if the global 
interaction is absent, three phases appear by increasing 
the interaction strength, the localized phase, delocalized 
chaotic phase and the phase of the ballistic motion. The 
final ballistic phase can be understood easily. When the 
direction of the velocity and hence the orientation of each 
particle is random, the average S is zero in eq. (fT2")l so 
that this equation is decoupled for each particle. The 
damping constant k is renormalized as n + K and hence, 
for large values of K, this is in the phase of straight 
motion in Fig. [H[a). It is expected qualitatively that 
the chaotic dynamics in the intermediate parameter re- 
gion originates from a competition between the circular 
motion of individual particles and the tendency of the 
orientation order for increasing the interaction strength. 
Finally it is mentioned that our preliminary numerical 
simulations have verified that the above property of the 
complex collective dynamics is preserved when an attrac- 
tive interaction whose strength is proportional to the dis- 
tance of a pair of particles. The full account of these re- 
sults together with the theoretical analysis will be pub- 
lished elsewhere in the neat future. 
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